/*! 
 ***************************************************************************
 * \file WaterFlooding1D.cpp
 * Two phase flow in porous media in 1D.
 *
 * \par Description
 * This code solves the 1D two phase flow in porous media in 1D, the
 * equations are posed in pressure-saturation formulation:
 * 
 * These equation are solved using the IMPES algorithm.
 * The domain of study is the interval \f$ x \in [0,0.25 m] \f$, and is
 * initially filled with oil. Then, water is injected on the left side
 * and we study the traveling front between oil and water.
 *
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./WaterFlooding1D  \endverbatim
 ***************************************************************************
 *  \author    Luis M. de la Cruz [ Apr 17 09:37:17 GMT 2011 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/TwoPhaseEquation.hpp"
#include "./WF/WFIP1.hpp"
#include "./WF/WFES1.hpp"
#include "Solvers/TDMA.hpp"
#include "Solvers/explicit.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,1>::huge ScalarField1D;

int pausa;

int main() 
{
    double length, dt, Tmax, tolerance; 
    int num_nodes, frequency;
    double porosity, permeability, mu_w, mu_o, Srw, Sro, theta, p_t;
    double injection, pres_out, Sw_in;

// Reading data from the input file
    std::ifstream input_cin ("inputWF1D"); 
    input_cin >> length
	      >> num_nodes
	      >> porosity
	      >> permeability
	      >> mu_w
	      >> mu_o
	      >> Srw
	      >> Sro
	      >> theta
	      >> p_t
	      >> injection
	      >> pres_out
	      >> Sw_in
	      >> dt
	      >> Tmax
	      >> frequency;
    input_cin.close();

    double Gamma = permeability;

    StructuredMesh<Uniform<double, 1> > mesh(length, num_nodes);
    float dx = mesh.getDelta(X);
    mesh.print();

    ScalarField1D p ( mesh.getExtentVolumes() );
    ScalarField1D Sw ( mesh.getExtentVolumes() ); 
    ScalarField1D p_n( mesh.getExtentNodes()   );
    ScalarField1D Sw_n( mesh.getExtentNodes()   );

    // Initial conditions
    Sw = Sw_in;  
    p = pres_out; 

    /* *
    std::cout << "\n Sw inicial = \n" << Sw;
    std::cout << "\n P inicial = \n" << p;
    std::cin >> pausa;
    /* */

    /* */
    NumUtils::interpolateToNodes(p_n, p);
    InOut::writeToFile(p_n, 0, "./DataWF1D/pres.", dx);
    NumUtils::interpolateToNodes(Sw_n, Sw);
    InOut::writeToFile(Sw_n, 0, "./DataWF1D/satu.", dx);

#ifdef WITH_GNUPLOT
    GNUplot plotter;
    int pause; 
    plotter("set grid");
#endif

    SparseMatrix< Diagonal< double, 1> > A(num_nodes); 
    ScalarField1D              b(num_nodes);

    TwoPhaseEquation< WFIP1<double, 1> > pressure(p, A, b, mesh.getDeltas());
    pressure.setDeltaTime(dt); 
    pressure.setPermeability(permeability);
    pressure.setPorosity(porosity);
    pressure.setSrw(Srw);
    pressure.setSro(Sro);
    pressure.setViscosity_w(mu_w);
    pressure.setViscosity_o(mu_o);
    pressure.setTheta(theta);
    pressure.setPthreshold(p_t);
    pressure.applyBounds(1, p.ubound(firstDim)-1); // just to move the indexes
    pressure.setPressure(p);
    pressure.setSaturation(Sw);
    pressure.print();

    //
    // In these case, the boundary conditions are calculated in
    // the numerical scheme (the adaptor BLES1).
    //
    TwoPhaseEquation< WFES1<double, 1> > saturation(Sw, A, b, mesh.getDeltas());
    saturation.setDeltaTime(dt); 
    saturation.setPermeability(permeability);
    saturation.setPorosity(porosity);
    saturation.setSrw(Srw);
    saturation.setSro(Sro);
    saturation.setViscosity_w(mu_w);
    saturation.setViscosity_o(mu_o);    
    saturation.setTheta(theta);
    saturation.setPthreshold(p_t);
    saturation.applyBounds(1, Sw.ubound(firstDim)-1); // just to move the indexes
    saturation.setPressure(p);
    saturation.setSaturation(Sw);
    saturation.print();

    /***/
    ofstream error_fp("./DataWF1D/error"), residual_fp("./DataWF1D/residual");
    int iteration = 0;
    double t = dt;
    double error = 1.0;
    double residual = 1.0;

    int ei = saturation.get_ei();

    while (t <= Tmax) {
      std::cout << "\n ---- Iteration " << ++iteration 
		<< "\t time = " << t
		<< "\t Res = " << residual
		<< "\t Err = " << error;
      
      pressure.calcCoefficients();

      Solver::TDMA1D(pressure);

      residual = pressure.calcResidual();
      error = pressure.calcErrorL1();
      //      error_fp << iteration << "\t" << error << std::endl;
      //      residual_fp << iteration << "\t" << residual << std::endl;	
       pressure.update();
     
      saturation.calcCoefficients();
      Solver::solExplicit(saturation);
      error = saturation.calcErrorL1();
      saturation.update();

      t += dt;
      if ( !(iteration % frequency) ) {
	NumUtils::interpolateToNodes(p_n, p);
	InOut::writeToFile(p_n, 1, "./DataWF1D/pres.",dx);
	NumUtils::interpolateToNodes(Sw_n, Sw);
	InOut::writeToFile(Sw_n, 1, "./DataWF1D/satu.",dx);

#ifdef WITH_GNUPLOT
	plotter("p [0:0.30][0:1] \"DataWF1D/satu.00001\" title \"Saturation\" w lp");
#endif
      }

    }

#ifdef WITH_GNUPLOT
    std::cout << "\n \n Press 1 and enter to see the Pressure ";
    std::cin >> pause;	
    plotter("p \"DataWF1D/pres.00001\" title \"Pressure\" w lp");
    std::cout << "\n \n Press 1 and enter to see the Saturation ";
    std::cin >> pause;	
    plotter("p [0:0.25][0:1] \"DataWF1D/satu.00001\" title \"Saturation\" w lp");
    std::cout << "\n \n Press 1 and enter to see finish ";
    std::cin >> pause;
#endif

    /** **/


    return 0;
}




